home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / msdos / viewers / pvquan16 / quant / heckbert.c < prev    next >
C/C++ Source or Header  |  1992-11-30  |  17KB  |  568 lines
  1. /************************************************************************
  2.  *                                                                      *
  3.  *                  Copyright (c) 1991, Frank van der Hulst             *
  4.  *                          All Rights Reserved                         *
  5.  *                                                                      *
  6.  * Authors:                                                             *
  7.  *        FvdH - Frank van der Hulst (Wellington, NZ)                   *
  8.  *                                                                      *
  9.  * Versions:                                                            *
  10.  *    V1.1 910626 FvdH - QUANT released for DBW_RENDER                  *
  11.  *    V1.2 911021 FvdH - QUANT released for PoV Ray                     *
  12.  *    V1.4 920303 FvdH - Ported to GNU                                  *
  13.  *    V1.6 921023 FvdH - Produce multi-image GIFs                       *
  14.  *                     - Port to OS/2 IBM C Set/2                       *
  15.  *                                                                      *
  16.  ************************************************************************/
  17. /*
  18.  * This software is copyrighted as noted below.  It may be freely copied,
  19.  * modified, and redistributed, provided that the copyright notice is
  20.  * preserved on all copies.
  21.  *
  22.  * There is no warranty or other guarantee of fitness for this software,
  23.  * it is provided solely "as is".  Bug reports or fixes may be sent
  24.  * to the author, who may or may not act on them as he desires.
  25.  *
  26.  * You may not include this software in a program or other software product
  27.  * without supplying the source, or without informing the end-user that the
  28.  * source is available for no extra charge.
  29.  *
  30.  * If you modify this software, you should include a notice giving the
  31.  * name of the person performing the modification, the date of modification,
  32.  * and the reason for such modification.
  33. */
  34. /*
  35.  * colorquant.c
  36.  *
  37.  * Perform variance-based color quantization on a "full color" image.
  38.  * Author:    Craig Kolb
  39.  *        Department of Mathematics
  40.  *        Yale University
  41.  *        kolb@yale.edu
  42.  * Date:    Tue Aug 22 1989
  43.  * Copyright (C) 1989 Craig E. Kolb
  44.  * $Id: colorquant.c,v 1.3 89/12/03 18:27:16 craig Exp $
  45.  *
  46.  * $Log:    colorquant.c,v $
  47.  *
  48.  * Revision 1.4  91/06/26  16:00:00  Frank van der Hulst
  49.  * Ported to Turbo C;
  50.  *       Call farmalloc rather than malloc
  51.  *       Virtual memory added to swap box data to/from disk
  52.  *       Rewritten in ANSI C
  53.  *       Removed call to QuantHistogram() from colorquant, to allow two
  54.  *       image files to create one palette
  55.  *       Changed QuantHistogram() to read from file, rather than from an
  56.  *       array
  57.  *          Changed format of palette to conform with the VGA palette
  58.  *
  59.  * Revision 1.3  89/12/03  18:27:16  craig
  60.  * Removed bogus integer casts in distance calculation in makenearest().
  61.  *
  62.  * Revision 1.2  89/12/03  18:13:12  craig
  63.  * FindCutpoint now returns FALSE if the given box cannot be cut.  This
  64.  * to avoid overflow problems in CutBox.
  65.  * "whichbox" in GreatestVariance() is now initialized to 0.
  66.  *
  67.  */
  68.  
  69. #ifdef __TURBOC__
  70. #include <mem.h>
  71. #define HUGE 1.79e308
  72. #endif
  73. #include <math.h>
  74.  
  75. #include "quant.h"
  76. #include "heckbert.h"
  77.  
  78. #define MAX(x,y)    ((x) > (y) ? (x) : (y))
  79.  
  80. static unsigned long        NPixels = 0L;            /* total # of pixels */
  81.  
  82. static int neighbours[MAXCOLORS];
  83.  
  84. /*
  85.  * Compute the histogram of the image as well as the projected frequency
  86.  * arrays for the first world-encompassing box.
  87.  * We compute both the histogram and the proj. frequencies of
  88.  * the first box at the same time to save a pass through the
  89.  * entire image.
  90.  * The projected frequency arrays of the largest box are zeroed out as
  91.  * as part of open_box_file(), called from main().
  92.  */
  93.  
  94. void QuantHistogram(Box *box)
  95. {
  96. unsigned long *rf, *gf, *bf;
  97. UCHAR pixel[3];
  98.  
  99.     rf = box->freq[RED];
  100.     gf = box->freq[GREEN];
  101.     bf = box->freq[BLUE];
  102.  
  103.     while (get_pixel(pixel)) {
  104.         rf[pixel[RED]]++;
  105.         gf[pixel[GREEN]]++;
  106.         bf[pixel[BLUE]]++;
  107.         Histogram[(((pixel[RED]<<INPUT_BITS)|pixel[GREEN])<<INPUT_BITS)|pixel[BLUE]]++;
  108.         NPixels++;
  109.     }
  110. }
  111.  
  112. /*
  113.  * Compute mean and weighted variance of the given box.
  114.  */
  115. void BoxStats(Box HUGE_PTR box)
  116. {
  117. int i, color;
  118. unsigned long *freq;
  119. double mean, var;
  120.  
  121.     if(box->weight == 0) {
  122.         box->weightedvar = 0;
  123.         return;
  124.     }
  125.  
  126.     box->weightedvar = 0.0;
  127.     for (color = 0; color < 3; color++) {
  128.         var = mean = 0;
  129.         i = box->low[color];
  130.         freq = &box->freq[color][i];
  131.         for (; i < box->high[color]; i++, freq++) {
  132.             mean += i * *freq;
  133.             var += i*i* *freq;
  134.         }
  135.         box->mean[color] = mean / box->weight;
  136.         box->weightedvar += var - box->mean[color]*box->mean[color]*box->weight;
  137.     }
  138.     box->weightedvar /= NPixels;
  139. }
  140.  
  141. /*
  142.  * Return the number of the box in 'boxes' with the greatest variance.
  143.  * Restrict the search to those boxes with indices between 0 and n-1.
  144.  */
  145. int GreatestVariance(int n)
  146. {
  147.     int i, whichbox = 0;
  148.     double max;
  149.     Box *box;
  150.  
  151.     max = -1;
  152.     for (i = 0; i < n; i++) {
  153.         box = get_box_tmp(i);
  154.         if (box->weightedvar > max) {
  155.             max = box->weightedvar;
  156.             whichbox = i;
  157.         }
  158.     }
  159.     return whichbox;
  160. }
  161.  
  162. /*
  163.  * Update projected frequency arrays for two boxes which used to be
  164.  * a single box.
  165.  */
  166. void UpdateFrequencies(Box HUGE_PTR box1, Box HUGE_PTR box2)
  167. {
  168. unsigned long myfreq, HUGE_PTR h;
  169. int b, g, r;
  170. int roff;
  171.  
  172.     memset(box1->freq[0], 0, IN_COLOURS * sizeof(unsigned long));
  173.     memset(box1->freq[1], 0, IN_COLOURS * sizeof(unsigned long));
  174.     memset(box1->freq[2], 0, IN_COLOURS * sizeof(unsigned long));
  175.  
  176.     for (r = box1->low[0]; r < box1->high[0]; r++) {
  177.         roff = r << INPUT_BITS;
  178.         for (g = box1->low[1];g < box1->high[1]; g++) {
  179.             b = box1->low[2];
  180.             h = Histogram + (((roff | g) << INPUT_BITS) | b);
  181.             for (; b < box1->high[2]; b++) {
  182.                 if ((myfreq = *h++) == 0)
  183.                     continue;
  184.                 box1->freq[0][r] += myfreq;
  185.                 box1->freq[1][g] += myfreq;
  186.                 box1->freq[2][b] += myfreq;
  187.                 box2->freq[0][r] -= myfreq;
  188.                 box2->freq[1][g] -= myfreq;
  189.                 box2->freq[2][b] -= myfreq;
  190.             }
  191.         }
  192.     }
  193. }
  194.  
  195. /*
  196.  * Compute the 'optimal' cutpoint for the given box along the axis
  197.  * indicated by 'color'.  Store the boxes which result from the cut
  198.  * in newbox1 and newbox2.
  199.  */
  200. int FindCutpoint(Box HUGE_PTR box, int color, Box HUGE_PTR newbox1, Box HUGE_PTR newbox2)
  201. {
  202.     double u, v, max;
  203.     int i, maxindex, minindex, cutpoint;
  204.     unsigned long optweight, curweight;
  205.  
  206.     if (box->low[color] + 1 == box->high[color])
  207.         return FALSE;    /* Cannot be cut. */
  208.     minindex = (int)((box->low[color] + box->mean[color]) * 0.5);
  209.     maxindex = (int)((box->mean[color] + box->high[color]) * 0.5);
  210.  
  211.     cutpoint = minindex;
  212.     optweight = box->weight;
  213.  
  214.     curweight = 0.;
  215.     for (i = box->low[color] ; i < minindex ; i++)
  216.         curweight += box->freq[color][i];
  217.     u = 0.;
  218.     max = -1;
  219.     for (i = minindex; i <= maxindex ; i++) {
  220.         curweight += box->freq[color][i];
  221.         if (curweight == box->weight)
  222.             break;
  223.         u += (double)(i * box->freq[color][i]) / box->weight;
  224.         v = ((double)curweight / (box->weight-curweight)) *
  225.                 (box->mean[color]-u)*(box->mean[color]-u);
  226.         if (v > max) {
  227.             max = v;
  228.             cutpoint = i;
  229.             optweight = curweight;
  230.         }
  231.     }
  232.     cutpoint++;
  233.     *newbox1 = *newbox2 = *box;
  234.     newbox1->weight = optweight;
  235.     newbox2->weight -= optweight;
  236.     newbox1->high[color] = cutpoint;
  237.     newbox2->low[color] = cutpoint;
  238.     UpdateFrequencies(newbox1, newbox2);
  239.     BoxStats(newbox1);
  240.     BoxStats(newbox2);
  241.  
  242.     return TRUE;    /* Found cutpoint. */
  243. }
  244.  
  245. /*
  246.  * Cut the given box.  Returns TRUE if the box could be cut, FALSE otherwise.
  247.  */
  248. int CutBox(Box HUGE_PTR box, Box HUGE_PTR newbox)
  249. {
  250.     int i;
  251.     double totalvar[3];
  252.     static Box newboxes[3][2];  /* Only used by CutBox, but don't want it on stack */
  253.  
  254.     if (box->weightedvar == 0. || box->weight == 0)
  255.         /*
  256.          * Can't cut this box.
  257.          */
  258.         return FALSE;
  259.  
  260.     /*
  261.      * Find 'optimal' cutpoint along each of the red, green and blue
  262.      * axes.  Sum the variances of the two boxes which would result
  263.      * by making each cut and store the resultant boxes for
  264.      * (possible) later use.
  265.      */
  266.     for (i = 0; i < 3; i++) {
  267.         if (FindCutpoint(box, i, &newboxes[i][0], &newboxes[i][1]))
  268.             totalvar[i] = newboxes[i][0].weightedvar +
  269.                 newboxes[i][1].weightedvar;
  270.         else
  271.             totalvar[i] = HUGE;
  272.     }
  273.  
  274.     /*
  275.      * Find which of the three cuts minimized the total variance
  276.      * and make that the 'real' cut.
  277.      */
  278.     if (totalvar[RED] <= totalvar[GREEN] &&
  279.         totalvar[RED] <= totalvar[BLUE]) {
  280.         *box = newboxes[RED][0];
  281.         *newbox = newboxes[RED][1];
  282.     } else if (totalvar[GREEN] <= totalvar[RED] &&
  283.          totalvar[GREEN] <= totalvar[BLUE]) {
  284.         *box = newboxes[GREEN][0];
  285.         *newbox = newboxes[GREEN][1];
  286.     } else {
  287.         *box = newboxes[BLUE][0];
  288.         *newbox = newboxes[BLUE][1];
  289.     }
  290.  
  291.     return TRUE;
  292. }
  293.  
  294. /*
  295.  * Iteratively cut the boxes.
  296.  */
  297. CutBoxes(int colors)
  298. {
  299. int curbox, varbox;
  300. Box *box = get_box(0);
  301.  
  302.     box->low[RED] = box->low[GREEN] = box->low[BLUE] = 0;
  303.     box->high[RED] = box->high[GREEN] = box->high[BLUE] = IN_COLOURS;
  304.     box->weight = NPixels;
  305.  
  306.     BoxStats(box);
  307.     free_box(0);
  308.  
  309.     printf("%d Boxes: cutting box: ", colors);
  310.     for (curbox = 1; curbox < colors; curbox++) {
  311.         printf("%3d\b\b\b", curbox);
  312.         varbox = GreatestVariance(curbox);
  313.         if (CutBox(get_box(varbox), get_box(curbox)) == FALSE)
  314.                 break;
  315.         free_box(curbox);
  316.         free_box(varbox);
  317.     }
  318.     printf("Done\n");
  319.  
  320.     return curbox;
  321. }
  322.  
  323. /*
  324.  * Make the centroid of "boxnum" serve as the representative for
  325.  * each color in the box.
  326.  */
  327. void SetRGBmap(int boxnum, Box *box, int bits)
  328. {
  329. int r, g, b;
  330.  
  331.     for (r = box->low[RED]; r < box->high[RED]; r++) {
  332.         for (g = box->low[GREEN]; g < box->high[GREEN]; g++) {
  333.             for (b = box->low[BLUE]; b < box->high[BLUE]; b++) {
  334.                 RGBmap[(((r<<bits)|g)<<bits)|b]=(char)boxnum;
  335.             }
  336.         }
  337.     }
  338. }
  339.  
  340. /*
  341.  * In order to minimize our search for 'best representative', we form the
  342.  * 'neighbors' array.  This array contains the number of the boxes whose
  343.  * centroids *might* be used as a representative for some color in the
  344.  * current box.  We need only consider those boxes whose centroids are closer
  345.  * to one or more of the current box's corners than is the centroid of the
  346.  * current box. 'Closeness' is measured by Euclidean distance.
  347.  */
  348. int getneighbors(int num, int colors)
  349. {
  350.     int i, j;
  351.     Box *bp;
  352.     double dist, LowR, LowG, LowB, HighR, HighG, HighB, ldiff, hdiff;
  353.  
  354.     bp = get_box_tmp(num);
  355.  
  356.     ldiff = bp->low[RED] - bp->mean[RED];
  357.     ldiff *= ldiff;
  358.     hdiff = bp->high[RED] - bp->mean[RED];
  359.     hdiff *= hdiff;
  360.     dist = MAX(ldiff, hdiff);
  361.  
  362.     ldiff = bp->low[GREEN] - bp->mean[GREEN];
  363.     ldiff *= ldiff;
  364.     hdiff = bp->high[GREEN] - bp->mean[GREEN];
  365.     hdiff *= hdiff;
  366.     dist += MAX(ldiff, hdiff);
  367.  
  368.     ldiff = bp->low[BLUE] - bp->mean[BLUE];
  369.     ldiff *= ldiff;
  370.     hdiff = bp->high[BLUE] - bp->mean[BLUE];
  371.     hdiff *= hdiff;
  372.     dist += MAX(ldiff, hdiff);
  373.  
  374.     dist = sqrt(dist);
  375.  
  376.     /*
  377.      * Loop over all colors in the colormap, the ith entry of which
  378.      * corresponds to the ith box.
  379.      *
  380.      * If the centroid of a box is as close to any corner of the
  381.      * current box as is the centroid of the current box, add that
  382.      * box to the list of "neighbors" of the current box.
  383.      */
  384.     HighR = (double)bp->high[RED] + dist;
  385.     HighG = (double)bp->high[GREEN] + dist;
  386.     HighB = (double)bp->high[BLUE] + dist;
  387.     LowR = (double)bp->low[RED] - dist;
  388.     LowG = (double)bp->low[GREEN] - dist;
  389.     LowB = (double)bp->low[BLUE] - dist;
  390.     for (i = j = 0; i < colors; i++) {
  391.         bp = get_box_tmp(i);
  392.         if (LowR <= bp->mean[RED] && HighR >= bp->mean[RED] &&
  393.             LowG <= bp->mean[GREEN] && HighG >= bp->mean[GREEN] &&
  394.             LowB <= bp->mean[BLUE] && HighB >= bp->mean[BLUE])
  395.             neighbours[j++] = i;
  396.     }
  397.  
  398.     return j;    /* Return the number of neighbors found. */
  399. }
  400.  
  401. /*
  402.  * Assign representative colors to every pixel in a given box through
  403.  * the construction of the NearestColor array.  For each color in the
  404.  * given box, we look at the list of neighbors passed to find the
  405.  * one whose centroid is closest to the current color.
  406.  */
  407. void makenearest(int boxnum, int nneighbors, int bits)
  408. {
  409.     int n, b, g, r;
  410.     double rdist, gdist, bdist, dist, mindist;
  411.     int which, *np;
  412.     Box *box, *bp;
  413.  
  414.     box = get_box(boxnum);
  415.  
  416.     for (r = box->low[RED]; r < box->high[RED]; r++) {
  417.         for (g = box->low[GREEN]; g < box->high[GREEN]; g++) {
  418.             for (b = box->low[BLUE]; b < box->high[BLUE]; b++) {
  419. /*
  420.  * The following two lines should be commented out if the RGBmap is going
  421.  * to be used for images other than the one given.
  422.  */
  423.                 if (Histogram[(((r<<bits)|g)<<bits)|b] == 0)
  424.                     continue;
  425.                 mindist = HUGE;
  426.                 /*
  427.                  * Find the colormap entry which is
  428.                  * closest to the current color.
  429.                  */
  430.                 np = neighbours;
  431.                 for (n = 0; n < nneighbors; n++, np++) {
  432.                     bp = get_box_tmp(*np);
  433.                     rdist = r - bp->mean[RED];
  434.                     gdist = g - bp->mean[GREEN];
  435.                     bdist = b - bp->mean[BLUE];
  436.                     dist = rdist*rdist + gdist*gdist + bdist*bdist;
  437.                     if (dist < mindist) {
  438.                         mindist = dist;
  439.                         which = *np;
  440.                     }
  441.                 }
  442.                 /*
  443.                  * The colormap entry closest to this
  444.                  * color is used as a representative.
  445.                  */
  446.                 RGBmap[(((r<<bits)|g)<<bits)|b] = which;
  447.             }
  448.         }
  449.     }
  450.     free_box(boxnum);
  451. }
  452. /*
  453.  * Form colormap and NearestColor arrays.
  454.  */
  455. void find_colors(int colors, int bits)
  456. {
  457. int i;
  458. int num;
  459.  
  460.     /*
  461.      * Form map of representative (nearest) colors.
  462.      */
  463.     printf("Mapping colours for box: ");
  464.     for (i = 0; i < colors; i++) {
  465.         printf("%3d\b\b\b", i);
  466.         /*
  467.          * Create list of candidate neighbors and
  468.          * find closest representative for each
  469.          * color in the box.
  470.          */
  471.         num = getneighbors(i, colors);
  472.         makenearest(i, num, bits);
  473.     }
  474.     printf("Done\n");
  475. }
  476.  
  477. /*
  478.  * Compute RGB to colormap index map.
  479.  */
  480. void ComputeRGBMap(int colors, int bits, int fast)
  481. {
  482. int i;
  483.  
  484.     if (fast) {
  485.         /*
  486.          * The centroid of each box serves as the representative
  487.          * for each color in the box.
  488.          */
  489.         for (i = 0; i < colors; i++)
  490.             SetRGBmap(i, get_box_tmp(i), bits);
  491.     } else
  492.         /*
  493.          * Find the 'nearest' representative for each pixel.
  494.          */
  495.         find_colors(colors, bits);
  496. }
  497.  
  498. /*
  499.  * Perform variance-based color quantization on a 24-bit image.
  500.  *
  501.  * Input consists of:
  502.  *    in_file    Pointer to file containing of red, green and blue pixel
  503.  *                intensities stored as unsigned characters.
  504.  *                The color of the ith pixel is given consecutive bytes, in the
  505.  *                order red, then green, then blue. Only the LS 4 bits are used.
  506.  *                0 indicates zero intensity, 15 full intensity.
  507.  *    pixels    The length of the red, green and blue arrays
  508.  *                in bytes, stored as an unsigned long int.
  509.  *    colormap    Points to the colormap.  The colormap
  510.  *                consists of red, green and blue arrays.
  511.  *                The red/green/blue values of the ith
  512.  *                colormap entry are given respectively by
  513.  *                colormap[0][i], colormap[1][i] and
  514.  *                colormap[2][i].  Each entry is an unsigned char.
  515.  *    colors    The number of colormap entries, stored
  516.  *                as an integer.
  517.  *    bits        The number of significant bits in each entry
  518.  *                of the red, green and blue arrays. An integer.
  519.  *    rgbmap    An array of unsigned chars of size (2^bits)^3.
  520.  *                This array is used to map from pixels to
  521.  *                colormap entries.  The 'prequantized' red,
  522.  *                green and blue components of a pixel
  523.  *                are used as an index into rgbmap to retrieve
  524.  *                the index which should be used into the colormap
  525.  *                to represent the pixel.  In short:
  526.  *                index = rgbmap[(((r << bits) | g) << bits) | b];
  527.  *     fast    If non-zero, the rgbmap will be constructed
  528.  *                quickly.  If zero, the rgbmap will be built
  529.  *                much slower, but more accurately.  In most
  530.  *                cases, fast should be non-zero, as the error
  531.  *                introduced by the approximation is usually
  532.  *                small.  'Fast' is stored as an integer.
  533.  *    Cfactor    Conversion factor.
  534.  *
  535.  * colorquant returns the number of colors to which the image was
  536.  * quantized.
  537.  */
  538.  
  539. int colorquant(int colors, int bits, int fast, double Cfactor)
  540. {
  541. int    i,                        /* Counter */
  542.         OutColors;            /* # of entries computed */
  543. Box    *box;
  544.  
  545.     OutColors = CutBoxes(colors);
  546.     /*
  547.      * We now know the set of representative colors.  We now
  548.      * must fill in the colormap and convert the representatives
  549.      * from their 'prequantized' range to 0-FULLINTENSITY.
  550.      */
  551.     for (i = 0; i < OutColors; i++) {
  552.         box = get_box_tmp(i);
  553.         palette[i][RED]   = (char)(box->mean[RED] * Cfactor + 0.5);
  554.         palette[i][GREEN] = (char)(box->mean[GREEN] * Cfactor + 0.5);
  555.         palette[i][BLUE]  = (char)(box->mean[BLUE] * Cfactor + 0.5);
  556.     }
  557.  
  558.     ComputeRGBMap(OutColors, bits, fast);
  559.  
  560.     return OutColors;        /* Return # of colormap entries */
  561. }
  562.  
  563. int pal_index(UCHAR *pixel)
  564. {
  565.     return RGBmap[(((pixel[RED]<<INPUT_BITS)|pixel[GREEN])<<INPUT_BITS)|pixel[BLUE]];
  566. }
  567.  
  568.